/*
 * fft_basic_redux.hpp
 *
 *  Created on: May 19, 2010
 *      Author: Jamie LaPointe
 */

#if !defined(FFT_BASIC_REDUX_HPP_)
#define FFT_BASIC_REDUX_HPP_

#include "fft_constants.hpp"
#include "psi_table.hpp"
#include <cmath>
#include <cstring>
#include <complex>

// <DEBUG>
#include <fstream>
#include <iostream>
// </DEBUG>

namespace FftLib
{
namespace BasicRedux
{

using namespace FftLib::Common;

class Fft
{
public:

   //
   // Generate a one dimensional bit reverse table
   //
   void generateBitReverse1dTable(int32 *table, int32 power)
   {
      int32 i = 0;
      int32 j = 0;
      int32 k = 0;
      int32 m = 0;
      int32 n = 0;
      int32 p = 0;

      table[0] = 0; // zeroth power
      for (i = 0; i < power; ++i)
      {
         n = 2 << i;
         p = n - 1;

         j = 0;
         for (k = 0; k < n; ++k)
         {
            table[p + k] = j;
            m = n / 2;

            while (m <= j && m != 0)
            {
               j -= m;
               m = m / 2;
            }

            j += m;
         }
      }
   }

   //
   // generate a two dimensional bit reverse table
   //
   void generateBitReverse2dTable(int32 *table2d, int32 *table1d, int32 m,
         int32 n)
   {
      int32 *tabr = table1d + m - 1;
      int32 *tabc = table1d + n - 1;
      int32 *tab2dr = table2d;
      int32 *tab2di = table2d + 1;
      int32 i = 0;
      int32 j = 0;
      int32 tn = 2 * n;
      int32 row = 0;
      int32 col = 0;
      int32 entryr = 0;
      int32 entryc = 0;

      for (i = 0; i < m; ++i)
      {
         row = i * tn;
         entryr = tabr[i] * tn;

         for (j = 0; j < n; ++j)
         {
            col = j * 2;
            entryc = tabc[j] * 2;

            tab2dr[row + col] = entryr + entryc;
            tab2di[row + col] = entryr + entryc;
         }
      }
   }

   //
   // generate a two dimensional Hanning filter
   //
   void generateHanning2dFilter(float32 *hanning, int32 m, int32 n)
   {
      float32 hrow[MaxHanningLength];
      float32 hcol[MaxHanningLength];
      float32 div = 0.0;
      float32 theta = 0.0;
      int32 i = 0;
      int32 j = 0;
      int32 r = 0;

      div = m + 1;
      for (i = 0; i < m; ++i)
      {
         theta = 2 * Pi * (i + 1) / div;
         hrow[i] = 0.5 * (1 - cosf(theta));
      }

      div = n + 1;
      for (i = 0; i < n; ++i)
      {
         theta = 2 * Pi * (i + 1) / div;
         hcol[i] = 0.5 * (1 - cosf(theta));
      }

      for (i = 0; i < m; ++i)
      {
         r = i * n;
         for (j = 0; j < n; ++j)
         {
            hanning[r + j] = hrow[i] * hcol[j];
         }
      }
   }

   //
   // Setup a two dimensional filter
   //
   void setup2dFilter(const float32 *kernel, int32 km, int32 kn,
         float32 *filter, int32 fm, int32 fn, float32 factor)
   {
      int32 i = 0;
      int32 j = 0;
      int32 ii = 0;
      int32 jj = 0;
      int32 tn = 2 * fn;
      int32 dir = 0;
      int32 djr = 0;
      int32 indm = 0;
      int32 indn = 0;

      memset(filter, 0, (2 * fm * fn) * sizeof(float32));

      dir = (km - 1) / 2;
      djr = (kn - 1) / 2;
      for (i = 0; i < km; ++i)
      {
         for (j = 0; j < kn; ++j)
         {
            indm = (fm + i - dir) % fm;
            indn = (fn + j - djr) % fn;

            ii = indm * tn;
            jj = indn * 2;
            filter[ii + jj] = factor * kernel[i * kn + j];
         }
      }
   }

   //
   // In place multiplication between a 2D complex array, stored linearly,
   // and itself.  The routine multiplies each data entry by the complex
   // conjugate of itself.  This effects a modulus squared of the
   // complex value, written back over itself.
   //
   void inPlaceMultiplication(float32 *data, int32 m, int32 n)
   {
      float32 *hr = data;
      float32 *hi = data + 1;
      float32 mr = 0.0;
      int32 i = 0;
      int32 j = 0;
      int32 k = 0;
      int32 tn = 2 * n;
      int32 tot = 2 * m * n;

      for (i = 0; i < tot; i += tn)
      {
         for (j = 0; j < tn; j += 2)
         {
            k = i + j;
            mr = (hr[k] * hr[k]) + (hi[k] * hi[k]);

            hr[k] = mr;
            hi[k] = 0.0;
         }
      }
   }

   //
   // In place multiplication between a 2D complex array, stored linearly, and
   // itself.  Data is two real 2D arrays, packed into the real and imaginary
   // components and then Fourier transformed.  Resultant Data is in bit-
   // reversed order.  In addition, a 2D filter transfer function is passed in
   // and multiplied to both Data components in the spectral domain.  the routine
   // extracts the individual Fourier components, then multiplies them point-
   // wise, conjuating one of them.
   //
   void inPlaceMultiplication(float32 *data, int32 *bitReverseTable, int32 m,
         int32 n, float32 *filter)
   {
      float32 *hr = data;
      float32 *hi = data + 1;
      float32 *fm2 = filter;
      float32 sr = 0.0;
      float32 si = 0.0;
      float32 ar = 0.0;
      float32 ai = 0.0;
      float32 pr = 0.0;
      float32 pi = 0.0;
      int32 *n1 = bitReverseTable;
      int32 i = 0;
      int32 j = 0;
      int32 k1 = 0;
      int32 k2 = 0;
      int32 tot = 2 * m * n;
      int32 htot = m * n;
      int32 tn = 2 * n;

      hr[0] = 4 * hr[0] * hi[0] * fm2[0];
      hi[0] = 0.0;

      // first row
      for (j = 2; j <= n; j += 2)
      {
         k1 = j;
         k2 = tn - j;
         sr = hr[n1[k1]] + hr[n1[k2]];
         si = hi[n1[k1]] - hi[n1[k2]];
         ar = hi[n1[k1]] + hi[n1[k2]];
         ai = hr[n1[k1]] - hr[n1[k2]];

         pr = sr * ar - si * ai;
         pi = sr * ai + si * ar;

         hr[n1[k1]] = pr * fm2[k1];
         hi[n1[k1]] = pi * fm2[k1];
         hr[n1[k2]] = pr * fm2[k2];
         hi[n1[k2]] = -pi * fm2[k2];
      }

      // first col
      for (i = tn; i <= m * n; i += tn)
      {
         k1 = i;
         k2 = tot - i;
         sr = hr[n1[k1]] + hr[n1[k2]];
         si = hi[n1[k1]] - hi[n1[k2]];
         ar = hi[n1[k1]] + hi[n1[k2]];
         ai = hr[n1[k1]] - hr[n1[k2]];

         pr = sr * ar - si * ai;
         pi = sr * ai + si * ar;

         hr[n1[k1]] = pr * fm2[k1];
         hi[n1[k1]] = pi * fm2[k1];
         hr[n1[k2]] = pr * fm2[k2];
         hi[n1[k2]] = -pi * fm2[k2];
      }

      // Main body
      for (i = tn; i < htot; i += tn)
      {
         for (j = 2; j < n; j += 2)
         {
            k1 = i + j;
            k2 = (tot - i) + (tn - j);
            sr = hr[n1[k1]] + hr[n1[k2]];
            si = hi[n1[k1]] - hi[n1[k2]];
            ar = hi[n1[k1]] + hi[n1[k2]];
            ai = hr[n1[k1]] - hr[n1[k2]];

            pr = sr * ar - si * ai;
            pi = sr * ai + si * ar;

            hr[n1[k1]] = pr * fm2[k1];
            hi[n1[k1]] = pi * fm2[k1];
            hr[n1[k2]] = pr * fm2[k2];
            hi[n1[k2]] = -pi * fm2[k2];
         }
      }
      for (i = htot; i < tot; i += tn)
      {
         for (j = 2; j <= n; j += 2)
         {
            k1 = i + j;
            k2 = (tot - i) + (tn - j);
            sr = hr[n1[k1]] + hr[n1[k2]];
            si = hi[n1[k1]] - hi[n1[k2]];
            ar = hi[n1[k1]] + hi[n1[k2]];
            ai = hr[n1[k1]] - hr[n1[k2]];

            pr = sr * ar - si * ai;
            pi = sr * ai + si * ar;

            hr[n1[k1]] = pr * fm2[k1];
            hi[n1[k1]] = pi * fm2[k1];
            hr[n1[k2]] = pr * fm2[k2];
            hi[n1[k2]] = -pi * fm2[k2];
         }
      }
   }

   //
   // Pack two real two dimensional arrays into one complex two dimensional
   // data array, stored linearly, with one real data array stored in the
   // real and the other stored in the imaginary parts.
   //
   void pack2dWindow(const int32 *in1ul, const int32 *in2ul, int32 cols,
         int32 subm, int32 subn, const float32 *window1,
         const float32 *window2, int32 padrad, float32 *packedData)
   {
      float32 *pdr = packedData;
      float32 *pdi = packedData + 1;
      int32 i = 0;
      int32 j = 0;
      int32 r = 0;
      int32 rp = 0;
      int32 cp = 0;
      int32 rr = 0;
      int32 tn = 2 * subn;

      // first, do the initial zero pad rows (2nd image)
      for (i = 0; i < padrad; ++i)
      {
         rp = i * tn;
         rr = i * subn;
         r = i * cols;
         for (j = 0; j < subn; ++j)
         {
            cp = j * 2;
            pdr[rp + cp] = window1[rr + j] * (float32) in1ul[r + j];
            pdi[rp + cp] = 0.0;
         }
      }

      // next, do the last zero pad rows
      for (i = subm - padrad; i < subm; ++i)
      {
         rp = i * tn;
         rr = i * subn;
         r = i * cols;
         for (j = 0; j < subn; ++j)
         {
            cp = j * 2;
            pdr[rp + cp] = window1[rr + j] * (float32) in1ul[r + j];
            pdi[rp + cp] = 0.0;
         }
      }

      // finally, do everything else
      for (i = padrad; i < subm - padrad; ++i)
      {
         rp = i * tn;
         rr = i * subn;
         r = i * cols;
         // do zero pad columns
         for (j = 0; j < padrad; ++j)
         {
            cp = j * 2;
            pdr[rp + cp] = window1[rr + j] * (float32) in1ul[r + j];
            pdi[rp + cp] = 0.0;
         }
         for (j = subn - padrad; j < subn; ++j)
         {
            cp = j * 2;
            pdr[rp + cp] = window1[rr + j] * (float32) in1ul[r + j];
            pdi[rp + cp] = 0.0;
         }

         // do the actual data
         for (j = padrad; j < subn - padrad; ++j)
         {
            cp = j * 2;
            pdr[rp + cp] = window1[rr + j] * (float32) in1ul[r + j];
            pdi[rp + cp] = window2[rr + j] * (float32) in2ul[r + j];
         }
      }
   }

   //
   // Decimation in frequency FFT using single precision with 'n' dimensions
   //
   // NOTE: this is NOT called every frame.
   void decimationInFrequency(float32 *data, const uint32 *nn, int32 ndim,
         int32 expsgn)
   {
      int32 idim = 0;
      uint32 i1 = 0;
      uint32 i2 = 0;
      uint32 i3 = 0;
      uint32 ip1 = 0;
      uint32 ip2 = 0;
      uint32 ip3 = 0;
      uint32 ifp1 = 0;
      uint32 ifp2 = 0;
      uint32 k1 = 0;
      uint32 k2 = 0;
      uint32 n = 0;
      uint32 nprev = 1;
      uint32 nrem = 0;
      uint32 ntot = 1;
      uint32 i2rev = 0;
      uint32 i3rev = 0;
      uint32 ibit = 0;
      float32 tr = 0.0;
      float32 trhn = 0.0;
      float32 ti = 0.0;
      float32 tihn = 0.0;
      float32 psi = 0.0;
      float32 wi = 0.0;
      float32 wni = 0.0;
      float32 wr = 0.0;
      float32 wnr = 0.0;
      float32 wt = 0.0;
      float32 *datar = data;
      float32 *datai = data + 1;

      for (idim = 0; idim < ndim; ++idim)
      {
         ntot = ntot * nn[idim];
      }

      for (idim = ndim - 1; idim >= 0; --idim)
      {
         n = nn[idim];
         nrem = ntot / (n * nprev);

         ip1 = nprev << 1;
         ip2 = ip1 * n;
         ip3 = ip2 * nrem;

         ifp1 = ip2;
         while (ifp1 > ip1)
         {
            ifp2 = ifp1 >> 1;
            psi = expsgn * (2 * Pi / (ifp1 / ip1));
            wnr = cosf(psi);
            wni = sinf(psi);

            wr = 1.0;
            wi = 0.0;
            for (i3 = 0; i3 < ifp2; i3 += ip1)
            {
               for (i2 = i3; i2 < i3 + ip1; i2 += 2)
               {
                  for (i1 = i2; i1 < ip3; i1 += ifp1)
                  {
                     k1 = i1;
                     k2 = k1 + ifp2;

                     tr = datar[k1];
                     ti = datai[k1];
                     trhn = datar[k2];
                     tihn = datai[k2];

                     datar[k1] = tr + trhn;
                     datai[k1] = ti + tihn;
                     datar[k2] = wr * (tr - trhn) - wi * (ti - tihn);
                     datai[k2] = wr * (ti - tihn) + wi * (tr - trhn);
                  }
               }
               wt = wr;
               wr = wr * wnr - wi * wni;
               wi = wt * wni + wi * wnr;
            }
            ifp1 = ifp2;
         }

         i2rev = 0;
         for (i2 = 0; i2 < ip2; i2 += ip1)
         {
            if (i2 < i2rev)
            {
               for (i1 = i2; i1 < i2 + ip1; i1 += 2)
               {
                  for (i3 = i1; i3 < ip3; i3 += ip2)
                  {
                     i3rev = i2rev + i3 - i2;
                     swap(&wt, &(datar[i3]), &(datar[i3rev]));
                     swap(&wt, &(datai[i3]), &(datai[i3rev]));
                  }
               }
            }
            ibit = ip2 >> 1;
            while (ibit <= i2rev)
            {
               i2rev -= ibit;
               ibit = ibit >> 1;
            }
            i2rev += ibit;
         }

         nprev *= n;
      }
   }

   //
   // Decimation in frequency FFT using single precision with 'n' dimensions
   // no bit reverse step is performed
   //
   // TODO: Investigate potential SIMD opportunities
   // TODO: Investigate other algorithm improvements based on known dimensions.
   //
   // This is called every frame when in CE
   //
   void decimationInFrequencyNoBitReverse(float32 *data, const uint32 *nn,
         int32 ndim, int32 expsgn)
   {
      int32 idim = 0;
      uint32 i1 = 0;
      uint32 i2 = 0;
      uint32 i3 = 0;
      uint32 ip1 = 0;
      uint32 ip2 = 0;
      uint32 ip3 = 0;
      uint32 ifp1 = 0;
      uint32 ifp2 = 0;
      uint32 k1 = 0;
      uint32 k2 = 0;
      uint32 n = 0;
      uint32 nprev = 1;
      uint32 nrem = 0;
      uint32 ntot = 1;
      float32 tr = 0.0;
      float32 trhn = 0.0;
      float32 ti = 0.0;
      float32 tihn = 0.0;
      float32 psi = 0.0;
      float32 wi = 0.0;
      float32 wni = 0.0;
      float32 wr = 0.0;
      float32 wnr = 0.0;
      float32 wt = 0.0;
      float32 *datar = data;
      float32 *datai = data + 1;

      for (idim = 0; idim < ndim; ++idim)
      {
         ntot = ntot * nn[idim];
      }

      for (idim = ndim - 1; idim >= 0; --idim)
      {
         n = nn[idim];
         nrem = ntot / (n * nprev);

         ip1 = nprev << 1;
         ip2 = ip1 * n;
         ip3 = ip2 * nrem;

         ifp1 = ip2;
         while (ifp1 > ip1)
         {
            ifp2 = ifp1 >> 1;
            // TODO: Investigate improvements to below calculations
            psi = expsgn * (2 * Pi / (ifp1 / ip1));
            wnr = cosf(psi);
            wni = sinf(psi);

            wr = 1.0;
            wi = 0.0;
            for (i3 = 0; i3 < ifp2; i3 += ip1)
            {
               for (i2 = i3; i2 < i3 + ip1; i2 += 2)
               {
                  for (i1 = i2; i1 < ip3; i1 += ifp1)
                  {
                     k1 = i1;
                     k2 = k1 + ifp2;

                     tr = datar[k1];
                     ti = datai[k1];
                     trhn = datar[k2];
                     tihn = datai[k2];

                     datar[k1] = tr + trhn;
                     datai[k1] = ti + tihn;
                     datar[k2] = wr * (tr - trhn) - wi * (ti - tihn);
                     datai[k2] = wr * (ti - tihn) + wi * (tr - trhn);
                  }
               }
               wt = wr;
               wr = wr * wnr - wi * wni;
               wi = wt * wni + wi * wnr;
            }
            ifp1 = ifp2;
         }

         nprev *= n;
      }
   }

   void create2dPsiTableDif(PsiTable *table, uint32 width, int32 expsgn)
   {
      int32 dimIndex = 0;
      uint32 ip1 = 0;
      uint32 ip2 = 0;
      uint32 ifp2 = 0;
      uint32 nprev = 1;
      uint32 index = 0.0;
      float32 psi = 0.0;

      for (dimIndex = 1; dimIndex >= 0; --dimIndex)
      {
         ip1 = nprev << 1; // 2, 128
         ip2 = ip1 * width; // 128, 8192

         // 64 > 2 (6 times), 8192 > 128 (6 times); [total of 12]
         while (ip2 > ip1)
         {
            ifp2 = ip2 >> 1;
            psi = expsgn * (TwoPi * ip1 / ip2);
            table[index].wnr = cosf(psi);
            table[index].wni = sinf(psi);
            ++index;
            ip2 = ifp2;
         }

         nprev *= width;
      }
   }

   // assumes a 2d square shaped matrix of data (width = width)
   void decimationInFrequencyNoBitReverse2d(float32 *data, uint32 width,
         PsiTable *psiTable)
   {
      int32 dimIndex = 0;
      uint32 i1 = 0;
      uint32 i2 = 0;
      uint32 i3 = 0;
      uint32 ip1 = 0;
      uint32 ip2 = 0;
      uint32 ip3 = 0;
      uint32 ifp1 = 0;
      uint32 ifp2 = 0;
      uint32 k1 = 0;
      uint32 k2 = 0;
      uint32 nprev = 1;
      uint32 nrem = 0;
      uint32 ntot = width * width; // 64x64 = 4096
      uint32 psiIndex = 0;
      float32 tr = 0.0;
      float32 trhn = 0.0;
      float32 ti = 0.0;
      float32 tihn = 0.0;
      float32 wi = 0.0;
      float32 wr = 0.0;
      float32 wt = 0.0;
      float32 *datar = data;
      float32 *datai = data + 1;

      for (dimIndex = 1; dimIndex >= 0; --dimIndex) // 1, 0
      {
         // nprev = 1, 64
         nrem = ntot / (width * nprev); // 64, 1

         ip1 = nprev << 1; //    2,  128
         ip2 = ip1 * width; //  128, 8192
         ip3 = ip2 * nrem; // 8192, 8192

         ifp1 = ip2; // 128, 8192
         // 128 > 2 (6 times), 8192 > 128 (6 times); total of 12 times
         while (ifp1 > ip1)
         {
            ifp2 = ifp1 >> 1;
            // TODO: Investigate improvements to below calculations

            wr = 1.0;
            wi = 0.0;
            for (i3 = 0; i3 < ifp2; i3 += ip1)
            {
               for (i2 = i3; i2 < i3 + ip1; i2 += 2)
               {
                  for (i1 = i2; i1 < ip3; i1 += ifp1)
                  {
                     k1 = i1;
                     k2 = k1 + ifp2;

                     tr = datar[k1];
                     ti = datai[k1];
                     trhn = datar[k2];
                     tihn = datai[k2];

                     datar[k1] = tr + trhn;
                     datai[k1] = ti + tihn;
                     datar[k2] = wr * (tr - trhn) - wi * (ti - tihn);
                     datai[k2] = wr * (ti - tihn) + wi * (tr - trhn);
                  }
               }
               wt = wr;
               wr = wr * psiTable[psiIndex].wnr - wi * psiTable[psiIndex].wni;
               wi = wt * psiTable[psiIndex].wni + wi * psiTable[psiIndex].wnr;
            }
            ++psiIndex;
            ifp1 = ifp2;
         }
         nprev *= width; // 1, 64
      }
   }

   // assume a 2x64x64 data matrix
   void decimationInFrequencyNoBitReverse4096(float32 *data, PsiTable *psiTable)
   {
      uint32 i1 = 0;
      uint32 i2 = 0;
      uint32 i3 = 0;
      uint32 ip2 = 0;
      uint32 ifp2 = 0;
      uint32 k1 = 0;
      uint32 k2 = 0;
      float32 tr = 0.0;
      float32 trhn = 0.0;
      float32 ti = 0.0;
      float32 tihn = 0.0;
      float32 wi = 0.0;
      float32 wr = 0.0;
      float32 wt = 0.0;
      float32 *datar = data;
      float32 *datai = data + 1;
      uint32 psiTableIndex = 0;

      ip2 = const4096::NumberOfRows; // 128

      while (ip2 > const4096::ComplexNumComponents) // 128 > 2 (6 times)
      {
         ifp2 = ip2 >> 1;
         // TODO: Investigate improvements to below calculations

         wr = 1.0;
         wi = 0.0;
         for (i3 = 0; i3 < ifp2; i3 += const4096::ComplexNumComponents)
         {
            for (i2 = i3; i2 < i3 + const4096::ComplexNumComponents; i2
                  += const4096::ComplexNumComponents)
            {
               for (i1 = i2; i1 < const4096::ComplexWindowArea; i1 += ip2)
               {
                  k1 = i1;
                  k2 = k1 + ifp2;

                  tr = datar[k1];
                  ti = datai[k1];
                  trhn = datar[k2];
                  tihn = datai[k2];

                  datar[k1] = tr + trhn;
                  datai[k1] = ti + tihn;
                  datar[k2] = wr * (tr - trhn) - wi * (ti - tihn);
                  datai[k2] = wr * (ti - tihn) + wi * (tr - trhn);
               }
            }
            wt = wr;
            wr = wr * psiTable[psiTableIndex].wnr - wi
                  * psiTable[psiTableIndex].wni;
            wi = wt * psiTable[psiTableIndex].wni + wi
                  * psiTable[psiTableIndex].wnr;
         }
         ++psiTableIndex;
         ip2 = ifp2;
      }

      ip2 = const4096::ComplexWindowArea; // 8192

      while (ip2 > const4096::NumberOfRows) // 8192 > 128 (6 times)
      {
         ifp2 = ip2 >> 1;
         // TODO: Investigate improvements to below calculations

         wr = 1.0;
         wi = 0.0;
         for (i3 = 0; i3 < ifp2; i3 += const4096::NumberOfRows)
         {
            for (i2 = i3; i2 < i3 + const4096::NumberOfRows; i2 += 2)
            {
               for (i1 = i2; i1 < const4096::ComplexWindowArea; i1 += ip2)
               {
                  k1 = i1;
                  k2 = k1 + ifp2;

                  tr = datar[k1];
                  ti = datai[k1];
                  trhn = datar[k2];
                  tihn = datai[k2];

                  datar[k1] = tr + trhn;
                  datai[k1] = ti + tihn;
                  datar[k2] = wr * (tr - trhn) - wi * (ti - tihn);
                  datai[k2] = wr * (ti - tihn) + wi * (tr - trhn);
               }
            }
            wt = wr;
            wr = wr * psiTable[psiTableIndex].wnr - wi
                  * psiTable[psiTableIndex].wni;
            wi = wt * psiTable[psiTableIndex].wni + wi
                  * psiTable[psiTableIndex].wnr;
         }
         ++psiTableIndex;
         ip2 = ifp2;
      }
   }

   //
   // Decimation in time FFT using single precision with 'n' dimensions no bit
   // reverse step is performed.
   //
   // TODO: Investigate potential SIMD opportunities
   // TODO: Investigate other algorithm improvements based on known dimensions.
   //
   // This is called every frame when in CE
   //
   void decimationInTimeNoBitReverse(float32 *data, const uint32 *nn,
         int32 ndim, int32 expsgn)
   {
      int32 idim = 0;
      uint32 i1 = 0;
      uint32 i2 = 0;
      uint32 i3 = 0;
      uint32 ip1 = 0;
      uint32 ip2 = 0;
      uint32 ip3 = 0;
      uint32 ifp1 = 0;
      uint32 ifp2 = 0;
      uint32 k1 = 0;
      uint32 k2 = 0;
      uint32 n = 0;
      uint32 nprev = 1; // 1, 64
      uint32 nrem = 0;
      uint32 ntot = 1;
      float32 tr = 0.0;
      float32 ti = 0.0;
      float32 psi = 0.0;
      float32 wi = 0.0;
      float32 wni = 0.0;
      float32 wnr = 0.0;
      float32 wr = 0.0;
      float32 wt = 0.0;
      float32 *datar = data;
      float32 *datai = data + 1;

      for (idim = 0; idim < ndim; ++idim)
      {
         ntot *= nn[idim]; // 4096
      }

      for (idim = ndim - 1; idim >= 0; --idim) // 1, 0
      {
         n = nn[idim]; // 64, 64
         nrem = ntot / (n * nprev); // 64, 1
         ip1 = nprev << 1; // 2, 128
         ip2 = ip1 * n; // 128, 8192
         ip3 = ip2 * nrem; // 8192, 8192

         ifp1 = ip1; // 2, 128
         while (ifp1 < ip2) // 2 < 128, 128 < 8192
         {
            ifp2 = ifp1 << 1;
            psi = expsgn * 2 * Pi / (ifp2 / ip1);
            wnr = cosf(psi);
            wni = sinf(psi);
            wr = 1.0;
            wi = 0.0;
            for (i3 = 0; i3 < ifp1; i3 += ip1)
            {
               for (i1 = i3; i1 < i3 + ip1; i1 += 2)
               {
                  for (i2 = i1; i2 < ip3; i2 += ifp2)
                  {
                     k1 = i2;
                     k2 = k1 + ifp1;
                     tr = wr * datar[k2] - wi * datai[k2];
                     ti = wr * datai[k2] + wi * datar[k2];
                     datar[k2] = datar[k1] - tr;
                     //datai[k2] = datai[k2] - ti; // I think this is WRONG!!!!!!!!!
                     datai[k2] = datai[k1] - ti;
                     datar[k1] += tr;
                     datai[k1] += ti;
                  }
               }
               wt = wr;
               wr = wr * wnr - wi * wni;
               wi = wi * wnr + wt * wni;
            }
            ifp1 = ifp2;
         }
         nprev *= n;
      }
   }

   void create2dPsiTableDit(PsiTable *table, uint32 width, int32 expsgn)
   {
      int32 dimIndex = 0;
      uint32 ip1 = 0;
      uint32 ifp1 = 0;
      uint32 ip2 = 0;
      uint32 ifp2 = 0;
      uint32 nprev = 1;
      uint32 index = 0.0;
      float32 psi = 0.0;

      for (dimIndex = 1; dimIndex >= 0; --dimIndex)
      {
         ip1 = nprev << 1; // 2, 128
         ip2 = ip1 * width; // 128, 8192
         ifp1 = ip1;

         // 2 < 128 (6 times), 128 < 8192 (6 times); [total of 12]
         while (ifp1 < ip2)
         {
            ifp2 = ifp1 << 1;
            psi = expsgn * (TwoPi * ip1 / ifp2);
            table[index].wnr = cosf(psi);
            table[index].wni = sinf(psi);
            ++index;
            ifp1 = ifp2;
         }

         nprev *= width;
      }
   }

   //
   // Decimation in time FFT using single precision with 'n' dimensions no bit
   // reverse step is performed.
   //
   // TODO: Investigate potential SIMD opportunities
   // TODO: Investigate other algorithm improvements based on known dimensions.
   //
   // This is called every frame when in CE
   //
   void decimationInTimeNoBitReverse4096(float32 *data, PsiTable *psiTable)
   {
      uint32 i1 = 0;
      uint32 i2 = 0;
      uint32 i3 = 0;
      uint32 ip1 = 0;
      uint32 ip2 = 0;
      uint32 ifp1 = 0;
      uint32 ifp2 = 0;
      uint32 k1 = 0;
      uint32 k2 = 0;
      uint32 psiIndex = 0;
      float32 ur = 0.0;
      float32 ui = 0.0;
      float32 vr = 0.0;
      float32 vi = 0.0;
      float32 wi = 0.0;
      float32 wr = 0.0;
      float32 wt = 0.0;
      float32 *datar = data;
      float32 *datai = data + 1;

      ip1 = 2;
      ip2 = const4096::NumberOfRows; // 128 (64 complex entries)

      ifp1 = ip1; // == 2
      while (ifp1 < ip2) // [2, 4, 8, 16, 32, 64] < 128 (6 times)
      {
         ifp2 = ifp1 << 1; // [4, 8, 16, 32, 64, 128]
         wr = 1.0;
         wi = 0.0;
         // 0 < [2, 4, 8, 16, 32, 64] step 2
         for (i3 = 0; i3 < ifp1; i3 += ip1)
         {
            // (0 < [2, 4, 8, 16, 32, 64] step 2) < 8192 step [4, 8, 16, 32, 64, 128]
            for (i2 = i3; i2 < const4096::ComplexWindowArea; i2 += ifp2)
            {
               k1 = i2;
               k2 = k1 + ifp1;

               ur = datar[k1];
               ui = datai[k1];

               vr = wr * datar[k2] - wi * datai[k2];
               vi = wr * datai[k2] + wi * datar[k2];

               datar[k1] = ur + vr;
               datai[k1] = ui + vi;

               datar[k2] = ur - vr;
               datai[k2] = ui - vi;
            }
            wt = wr;
            wr = wr * psiTable[psiIndex].wnr - wi * psiTable[psiIndex].wni;
            wi = wi * psiTable[psiIndex].wnr + wt * psiTable[psiIndex].wni;
         }
         ++psiIndex;
         ifp1 = ifp2;
      }

      // **********************************************
      // ******* 2nd go around ************************
      // **********************************************

      ip1 = const4096::NumberOfRows; // 128 (64 complex items)
      ip2 = const4096::ComplexWindowArea; // 8192 (4096 complex items)

      ifp1 = ip1; // 128
      // [128, 256, 512, 1024, 2048, 4096] < 8192 (6 times)
      while (ifp1 < ip2)
      {
         ifp2 = ifp1 << 1; // [256, 512, 1024, 2048, 4096, 8192]
         wr = 1.0;
         wi = 0.0;
         // 0 < [128, 256, 512, 1024, 2048, 4096] step 128
         for (i3 = 0; i3 < ifp1; i3 += ip1)
         {
            // i3 < (i3 + 128) step 2 (64 times)
            for (i1 = i3; i1 < i3 + ip1; i1 += 2)
            {
               // i1 < 8129 step [256, 512, 1024, 2048, 4096, 8192]
               for (i2 = i1; i2 < const4096::ComplexWindowArea; i2 += ifp2)
               {
                  k1 = i2;
                  k2 = k1 + ifp1;

                  ur = datar[k1];
                  ui = datai[k1];

                  vr = wr * datar[k2] - wi * datai[k2];
                  vi = wr * datai[k2] + wi * datar[k2];

                  datar[k1] = ur + vr;
                  datai[k1] = ui + vi;

                  datar[k2] = ur - vr;
                  datai[k2] = ui - vi;
               }
            }
            wt = wr;
            wr = wr * psiTable[psiIndex].wnr - wi * psiTable[psiIndex].wni;
            wi = wi * psiTable[psiIndex].wnr + wt * psiTable[psiIndex].wni;
         }
         ++psiIndex;
         ifp1 = ifp2;
      }
   }

   //
   // Decimation in time FFT using single precision with 'n' dimensions no bit
   // reverse step is performed.
   //
   // TODO: Investigate potential SIMD opportunities
   // TODO: Investigate other algorithm improvements based on known dimensions.
   //
   // This is called every frame when in CE
   //
   void decimationInTimeNoBitReverse4096Complex(std::complex<float32> *data,
         PsiTable *psiTable)
   {
      uint32 i1 = 0;
      uint32 i2 = 0;
      uint32 i3 = 0;
      uint32 ip1 = 0;
      uint32 ip2 = 0;
      uint32 ifp1 = 0;
      uint32 ifp2 = 0;
      uint32 k1 = 0;
      uint32 k2 = 0;
      uint32 psiIndex = 0;
      std::complex<float32> u;
      std::complex<float32> v;
      std::complex<float32> w;
      std::ofstream myfile;

      ip1 = 1;
      ip2 = const4096::WindowWidth; // 64 complex entries

      ifp1 = ip1; // == 1
      while (ifp1 < ip2) // [1, 2, 4, 8, 16, 32] < 64 (6 times)
      {
         ifp2 = ifp1 << 1; // [2, 4, 8, 16, 32, 64]
         w = std::complex<float32>(1.0, 0.0);
         // 0 < [1, 2, 4, 8, 16, 32] step 1
         for (i3 = 0; i3 < ifp1; i3 += ip1)
         {
            // (0 < 1, 2, 4, 8, 16, 32] step 1) < 4096 step [2, 4, 8, 16, 32, 64]
            for (i2 = i3; i2 < const4096::WindowArea; i2 += ifp2)
            {
               k1 = i2;
               k2 = k1 + ifp1;

               u = data[k1];
               v = data[k2] * w;

               data[k1] = u + v;
               data[k2] = u - v;
            }
            std::complex<float32> psiTableComplex = std::complex<float32>(
                  psiTable[psiIndex].wnr, psiTable[psiIndex].wni);
            w = w * psiTableComplex;
         }
         ++psiIndex;
         ifp1 = ifp2;
      }

      myfile << std::endl << std::endl;

      // **********************************************
      // ******* 2nd go around ************************
      // **********************************************

      ip1 = const4096::WindowWidth; // 64 complex items
      ip2 = const4096::WindowArea; // 4096 complex items

      ifp1 = ip1; // 64
      // [64, 128, 256, 512, 1024, 2048] < 4096 (6 times)
      while (ifp1 < ip2)
      {
         ifp2 = ifp1 << 1; // [128, 256, 512, 1024, 2048, 4096]
         //wr = 1.0;
         //wi = 0.0;
         w = std::complex<float32>(1.0, 0.0);
         for (i3 = 0; i3 < ifp1; i3 += ip1)
         {
            // i3 < (i3 + 64) step 2 (64 times)
            for (i1 = i3; i1 < i3 + ip1; i1 += 1)
            {
               // i1 < 4096 step [128, 256, 512, 1024, 2048, 4096]
               for (i2 = i1; i2 < const4096::WindowArea; i2 += ifp2)
               {
                  k1 = i2;
                  k2 = k1 + ifp1;

                  u = data[k1];
                  v = data[k2] * w;

                  data[k1] = u + v;
                  data[k2] = u - v;
               }
            }
            std::complex<float32> psiTableComplex = std::complex<float32>(
                  psiTable[psiIndex].wnr, psiTable[psiIndex].wni);
            w = w * psiTableComplex;
         }
         ++psiIndex;
         ifp1 = ifp2;
      }
   }

private:
   static const int32 MaxHanningLength = 128;

   inline void swap(float32 *t, float32 *a, float32 *b)
   {
      *t = *a;
      *a = *b;
      *b = *t;
   }
};

}
}

#endif /* FFT_BASIC_REDUX_HPP_ */
